Zmiany złożoności przestrzennej obszarów leśnych w latach 1992-2020

Author

Adrian Nowacki

Wprowadzenie

Poniższy raport obejmuje przygotowanie danych oraz przebieg analizy i wyniki przygotowane dla jednej z trzech klas pokrycia terenu - obszarów leśnych. Posłużono się jedną, przykładową klasą, ze względu na skupieniu się na przedstawieniu metodyki oraz prezentacji głównych wyników. Celem analizy jest określenie zmian złożoności przestrzennej obszarów leśnych na świecie w latach 1992-2020.

Strukturę przestrzenną lokalnego krajobrazu zaprezentowanego w formie kategorycznego rastra można określić m. in. za pomocą macierzy współwystępowania, czyli relacji pomiędzy sąsiadującymi komórkami rastra. Na ich podstawie można obliczyć dwie główne metryki: entropię brzegową oraz względną informację wzajemną. Określają one kolejno różnorodność (złożoność) oraz konfigurację (jednolitość, zbitość) lokalnego krajobrazu. Czym większa wartość entropii brzegowej, tym większa różnorodność krajobrazu, czym większa wartość względnej informacji wzajemnej - tym krajobraz jest bardziej jednolity przestrzennie.


Przygotowanie danych

1.1 Wczytanie danych rastrowych

Rastry obejmujące globalne pokrycie terenu na lata 1992-2020 pozyskano z projektu ESA Climate Change Initiative (1992-2015) oraz Copernicus Climate Change Service (2016-2020). Charakteryzują się one rozdzielczością przestrzenną równą 300x300 metrów. Ze względu na bardzo znikome lub zerowe zmiany pokrycia terenu, Antarktyda została wyłączona z analizy.

folders = c("africa", "asia", "europe", "north_america", "oceania", "south_america")
esa_files <- list(paste0("ESACCI-LC-L4-LCCS-Map-300m-P1Y-", 1992:2020, "-v2.0.7"))

for (i in 1:length(esa_files)) {
  file_name <- paste0("esa_cci/", esa_files[i], ".tif")
  rast <- rast(file_name)
  year <- regmatches(esa_files[[i]], regexpr("\\d{4}", esa_files[[i]]))
  assign(paste0("LC_", year), rast)
}

1.2 Wczytanie danych wektorowych

Aby w pełni poprawnie obliczyć zmiany powierzchniowe, posłużono się danymi pochodzącymi z Equi7Grid - projektu odpowiadającego za utworzenie wiernoodległościowych systemów odniesienia dla wszystkich 7 kontynentów. Indywidualne siatki odwzorowań mają na celu zwiększenie efektywności przetwarzania wysokorozdzielczych danych przestrzennych.

for (i in 1:6) {
  files = list.files(paste0(folders[i]), pattern = "\\GEOG_LAND.shp$", full.names = TRUE)
  continents = vect(files)
  continents_proj = project(continents, "+proj=igh +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
  assign(folders[i], continents)
}

cont_layers = list(africa, asia, europe, north_america, oceania, south_america)
layers_crs <- c("+proj=aeqd +lat_0=8.5 +lon_0=21.5 +x_0=5621452.01998 +y_0=5990638.42298 +datum=WGS84 +units=m +no_defs", # africa
                "+proj=aeqd +lat_0=47 +lon_0=94 +x_0=4340913.84808 +y_0=4812712.92347 +datum=WGS84 +units=m +no_defs", # asia
                "+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +datum=WGS84 +units=m +no_defs", # europe
                "+proj=aeqd +lat_0=52 +lon_0=-97.5 +x_0=8264722.17686 +y_0=4867518.35323 +datum=WGS84 +units=m +no_defs", # north america
                "+proj=aeqd +lat_0=-19.5 +lon_0=131.5 +x_0=6988408.5356 +y_0=7654884.53733 +datum=WGS84 +units=m +no_defs", # oceania
                "+proj=aeqd +lat_0=-14 +lon_0=-60.5 +x_0=7257179.23559 +y_0=5592024.44605 +datum=WGS84 +units=m +no_defs") # south america

year <- 1992:2020
LC_list <- list(paste0('LC_', 1992:2020))

2. Przycięcie oraz przekształcenie danych rastrowych

Założeniem przygotowania danych było uzyskanie warstw rastrowych pokrycia terenu każdego kontynentu dla każdego roku. Głównym etapem przygotowania danych do analizy było zatem przycięcie rastrów do obszarów poszczególnych kontynentów oraz zmiana układów współrzędnych z WGS84 na indywidualne odwzorowania dostarczone przez Equi7Grid.

Z powodu długotrwałego przekształcenia danych przy użyciu pakietów raster oraz terra, zdecydowano się wykorzystać w tym celu gdalUtils, który znacznie zmniejszył czas obliczeń. Etap ten był najdłużej trwającym etapem przygotowywania danych i zajął łącznie kilkadziesiąt godzin przetwarzania, głównie ze względu na problematyczny system odwzorowania dla Oceanii. W celu usprawnienia obliczeń zdecydowano się na podzielenie rastra Oceanii na dwie części, według równoleżnika 180°, co znacznie zmniejszyło czas transformacji układu z ponad 5 godzin do ok. 40 minut dla tego kontynentu.

W pracy magisterskiej zapewne warto umieścić tabelę przedstawiającą czas, jaki zajęły poszczególne procesy.

transf_oceania <- function(){
      ## ------------ first part -----------
      rast_1 = crop(x, extent(70, 180, -60, 30))
      raster_cropped_1 <- paste0("projected/", folders[j], "_temp1", year, ".tif")
      writeRaster(rast_1, filename = raster_cropped_1, datatype = "INT1U")
      
      gdalwarp(srcfile = raster_cropped_1,
               dstfile = paste0("projected/", folders[j], "_proj0", year, ".tif"),
               t_srs = layers_crs[[j]], tr = c(300, 300), r = "near", wm = 1000)
      
      gdal_translate(paste0("projected/", folders[j], "_proj0", year, ".tif"), paste0("projected/", folders[j], "_proj_1_", year, ".tif"),
                     of = "COG", co = "compress=lzw")
      
      
      ## ---------- second part ---------
      rast_2 = crop(x, extent(-180, -125, -50, 10))
      raster_cropped_2 <- paste0("projected/", folders[j], "_temp2", year, ".tif")
      writeRaster(rast_2, filename = raster_cropped_2, datatype = "INT1U")
      gdalwarp(srcfile = raster_cropped_2,
               dstfile = paste0("projected/", folders[j], "_proj1", year, ".tif"),
               t_srs = layers_crs[[j]], tr = c(300, 300), r = "near", wm = 1000)
      gdal_translate(paste0("projected/", folders[j], "_proj1", year, ".tif"), paste0("projected/", folders[j], "_proj_2_", year, ".tif"),
                     of = "COG", co = "compress=lzw")
}


transf_cont <- function(){
      raster_cropped <- paste0("projected/", folders[j], "_temp.tif")
      writeRaster(x, filename = raster_cropped, datatype = "INT1U")
      gdalwarp(srcfile = raster_cropped,
               dstfile = paste0("projected/", folders[j], "_proj0", year, ".tif"),
               t_srs = layers_crs[[j]], tr = c(300, 300), r = "near", wm = 1000)
      gdal_translate(paste0("projected/", folders[j], "_proj0", year, ".tif"), paste0("projected/", folders[j], "_proj_", year, ".tif"),
                     of = "COG", co = "compress=lzw")
}

Przygotowano dwie osobne funkcje transf_oceania oraz transf_cont, odpowiadające za zmianę układu współrzędnych przyciętych rastrów. Funkcje te wykorzystano w poniższej pętli dla wszystkich plików rastrowych oraz wszystkich lat.

for(i in LC_list) {
  for (j in 1:length(cont_layers)) {
    x = crop(i, cont_layers[j])
    x = terra::mask(x, vect(cont_layers[j]))
    
    if(j == 5){
      transform_oceania()
      
      file.remove(paste0("projected/", folders[j], "_proj0", year, ".tif"))
      file.remove(paste0("projected/", folders[j], "_proj1", year, ".tif"))
      file.remove(raster_cropped_1)
      file.remove(raster_cropped_2)
      
      oceania_1 <- rast(paste0("projected/", folders[j], "_proj_1_", year, ".tif"))
      oceania_2 <- rast(paste0("projected/", folders[j], "_proj_2_", year, ".tif"))
      merged <- vrt(c(oceania_1, oceania_2), paste0("projected/", folders[j], "_proj_", year, ".tif"))
      
      file.remove(paste0("projected/", folders[j], "_proj_1_", year, ".tif"))
      file.remove(paste0("projected/", folders[j], "_proj_2_", year, ".tif"))
    } else{
      
      transf_cont()
      
      file.remove(paste0("projected/", folders[j], "_proj0", year, ".tif"))
      file.remove(raster_cropped)
    }
  }
  year <<- year + 1
}

3. Reklasyfikacja

Obliczenie zmian złożoności przestrzennej zawartej w temacie pracy magisterskiej wymagało trzykrotnej reklasyfikacji danych na dwie klasy: obszary leśne oraz pozostałe, obszary zurbanizowane oraz pozostałe, a także obszary rolnicze oraz pozostałe. Zabieg ten pozwolił na szczegółowe zbadanie zmian struktury przestrzennej lokalnych krajobrazów oraz głównych trendów, m. in.: wylesiania, ekspansji rolnictwa oraz urbanizacji.

Reklasyfikację przeprowadzono w oparciu o zgrupowane klasy pokrycia terenu, proponowane przez IPCC (str. 30) w celu uniwersalnego porównania zmian 6 głównych klas pokrycia terenu, a także w oparciu o dokumentację klas ESA CCI.

m_forest <- c(0, 49, 0,
              49, 101, 1,
              101, 159, 0,
              159, 171, 1,
              171, 221, 0
)
matrix_forest <- matrix(m_forest, ncol = 3, byrow = TRUE)

years <- 1992:2020
folders = c("africa", "asia", "europe", "north_america", "oceania", "south_america")

for (i in 1:length(years)) {
  files <- list.files(paste0(years[i]), pattern = ".tif$", full.names = TRUE)
  
  for (j in 1:length(files)) {
    continent <- rast(files[j])
    raster_forest <- classify(continent, matrix_forest, include.lowest = TRUE, brackets=FALSE)
    writeRaster(raster_forest, paste0("reclassified/", folders[j], "_forest_", years[i], ".tif"), datatype = "INT1U", overwrite = TRUE)
  }
}

Dane otrzymane po reklasyfikacji są zbiorem rastrów w podziale na lata, kontynenty oraz klasy pokrycia terenu.

library(terra)
afr_1992 <- rast("reclassified/1992/south_america_forest_1992.tif")
afr_2020 <- rast("reclassified/2020/south_america_forest_2020.tif")

plot(c(afr_1992, afr_2020), main = c("1992", "2020"), col = c("#444444", "#468b69"))

Przykładowe dane po reklasyfikacji: obszary leśne Ameryki Południowej w 1992 i 2020 roku


Przebieg analizy

Dalsze etapy pracy obejmują wykorzystanie zreklasyfikowanych danych do obliczenia metryk złożoności przestrzennej oraz udziału procentowego poszczególnych klas pokrycia terenu, w podziale na kontynenty oraz lata.

Obliczenie metryk

Wartości entropii brzegowej oraz względnej informacji wzajemnej otrzymano za pomocą kluczowych funkcji z pakietu motif. W celu ujednolicenia wyników, a także metod obliczenia statystyk i wykonania map, utworzono obszerną ramkę danych zawierającą niezbędne wyniki oraz informacje. Metryki wyliczono dla ponad 160 tys. lokalnych krajobrazów dla całego świata, każdy o rozdzielczości przestrzennej 30x30 km (100 razy mniejszej od rozdzielczości danych surowych). Otrzymana ramka metrics zawiera informacje o ID lokalnego krajobrazu, wartości entropii brzegowej, informacji wzajemnej oraz względnej informacji wzajemnej, a także o przyjętej rozdzielczości przestrzennej oraz procentowym udziale wartości NA, klasy pokrycia terenu i obszarów pozostałych. Dla każdego z nich dołączono również jego geometrię, uzyskując obszerny obiekt sf.

calc_relmutinf = function(mutinf, ent){
  ifelse(mutinf == 0, 1, mutinf / ent)
}

add_it_metrics = function(x, year, res){
  x %>%
    mutate(ent = map_dbl(signature, comat:::rcpp_ent)) %>% 
    mutate(mutinf = map_dbl(signature, comat:::rcpp_mutinf)) %>%
    mutate(relmutinf = calc_relmutinf(mutinf, ent)) %>% 
    mutate(res = res) %>%
    select(-signature) %>% 
    set_names(c("id", "naprop", "ent", "mutinf", "relmutinf", "res"))
}

calculate_all_metrics = function(year, continent, class){
  input_data = read_stars(file)
  window = 100
  res = 30
  
  input_signature = lsp_signature(input_data, type = "coma", window = window)
  input_metrics = add_it_metrics(input_signature, year, res = res)
  input_metrics <- lsp_add_sf(input_metrics)
  
  input_proc = lsp_signature(input_data, type = "composition", window = window)
  input_proc = lsp_restructure(input_proc)
  input_proc_sf = lsp_add_sf(input_proc)
  metrics <- cbind(input_metrics, 
                   "prop_0" = input_proc$X1,
                   "prop_1" = input_proc$X2,
                   "year" = year,
                   "continent" = continent,
                   "class" = class)
  names(metrics)[names(metrics) == "naprop"] <- "prop_NA"
  names(metrics)[names(metrics) == "prop_0"] <- "prop_0"
  names(metrics)[names(metrics) == "prop_1"] <- "prop_1"
  metrics
}

W celu jak najdokładniejszej identyfikacji lokalnych krajobrazów na późniejszych etapach analizy, na podstawie nazewnictwa rastrów dodano informację o kontynencie, roku oraz reprezentowanej klasie pokrycia terenu. Poniższą pętlę powtórzono dla każdej klasy pokrycia terenu, co łącznie przełożyło się na ponad 35 godzin przetwarzania danych.

years <- 1992:2020
metrics <- data.frame()

for (i in 1:length(years)) {
  files <- list.files(paste0(years[i]), pattern = "_forest_.*\\.tif$", full.names = TRUE)
  
  for (file in files) {
    file_name <- basename(file)
    num_tokens <- length(unlist(strsplit(file_name, "_")))
    
    if (num_tokens == 4) {
      continent  <- gsub("^([^_]+_[^_]+)_.*", "\\1", file_name)
      class <- gsub("^.*?_[^_]+_([^_]+)_.*", "\\1", file_name)
    } 
    else if (num_tokens == 3) {
      continent  <- gsub("^(.*?)_.*", "\\1", file_name)
      class <- gsub("^.*?_(.*?)_.*", "\\1", file_name)
    }
    year <- as.integer(gsub("^.*?_.*_(\\d{4}).*", "\\1", file_name))
    df <- calculate_all_metrics(year, continent, class)
    df <- as.data.frame(df)
    metrics <- rbind(metrics, df)
  }
}
metrics <- sample_n(metrics, 10)
subset(metrics, select = -geometry)
id prop_NA ent mutinf relmutinf res prop_0 prop_1 year continent class
23038 0.009 0.0000 0.0000 1.0000 30 1.0000 0.0000 1992 africa forest
80327 0.000 0.0000 0.0000 1.0000 30 1.0000 0.0000 2013 africa urban
60486 0.000 0.3325 0.2172 0.6531 30 0.9393 0.0607 2001 africa agro
53615 0.000 0.2559 0.0851 0.3325 30 0.0432 0.9568 1993 africa forest
439842 0.000 0.6787 0.4010 0.5908 30 0.8204 0.1796 2011 asia forest
331751 0.000 0.1012 0.0380 0.3751 30 0.9868 0.0132 2007 asia agro
268588 0.000 0.0000 0.0000 1.0000 30 1.0000 0.0000 2008 north_america agro
359020 0.000 0.0000 0.0000 1.0000 30 1.0000 0.0000 2000 asia urban
369952 0.000 0.0000 0.0000 1.0000 30 1.0000 0.0000 2018 asia urban
125687 0.000 0.0291 0.0162 0.5586 30 0.9969 0.0031 2018 oceania urban

Grupowanie metryk

W celu obliczenia zmian wartości metryk, z otrzymanego obiektu sf metrics wydzielono lata 1992 i 2020 oraz pogrupowano według kontynentu, klasy i geometrii, aby dodać kolumnę z unikalnym ID każdego lokalnego krajobrazu. Pogrupowane rekordy dla różnych lat przyjęły identyczne ID - to pozwoliło na obliczenie zmiany wartości entropii brzegowej oraz względnej informacji wzajemnej na podstawie par krajobrazów lokalnych. W celu ich późniejszej identyfikacji dodano informację o kontynencie i klasie pokrycia terenu, a także pierwotnym ID krajobrazu.

sf_result <- metrics %>%
  group_by(continent, class, geometry) %>%
  mutate(group_id = cur_group_id())

df_grouped <- sf_result %>%
  group_by(group_id) %>%
  summarize(ent_diff = diff(ent),
            mutinf_diff = diff(mutinf),
            relmutinf_diff = diff(relmutinf))

df_grouped$continent <- sf_result$continent[match(df_grouped$group_id, sf_result$group_id)]
df_grouped$class <- sf_result$class[match(df_grouped$group_id, sf_result$group_id)]
df_grouped$id <- sf_result$id[match(df_grouped$group_id, sf_result$group_id)]
df_grouped <- st_drop_geometry(df_grouped)
df_grouped <- df_grouped[order(df_grouped$id),]
df_grouped <- df_grouped[1:9,]
group_id ent_diff mutinf_diff relmutinf_diff continent class id
419470 -0.0473 -0.0103 -0.0147 south_america agro 175
440163 0.0691 0.0545 0.0465 south_america forest 175
460856 0.1993 0.0433 -0.0690 south_america urban 175
419471 0.0126 0.0008 -0.0043 south_america agro 788
440164 -0.0002 -0.0074 -0.0073 south_america forest 788
460857 0.0162 0.0090 0.4089 south_america urban 788
419472 0.0425 0.0619 0.0625 south_america agro 1088
440165 0.0224 0.0693 0.0641 south_america forest 1088
460858 0.0852 0.0380 0.0846 south_america urban 1088

Obliczenie powierzchni klas

W celu otrzymania szczegółowych zmian powierzchni, zreklasyfikowane rastry zamieniono na ramki danych, dodając do każdego piksela informacje identyfikacyjne. Ramkę tę następnie pogrupowano według kontynentu i roku, aby zsumować liczbę występowań pikseli i obliczyć udział procentowy każdej klasy.

years <- 1992:2020
result_df <- data.frame()

for (i in 1:length(years)) {
  files <- list.files(paste0(years[i]), pattern = ".tif$", full.names = TRUE)
  
  for (file in files) {
    raster_data <- rast(file)
    freq_data <- as.data.frame(freq(raster_data))
    file_name <- basename(file)
    
    num_tokens <- length(unlist(strsplit(file_name, "_")))
    if (num_tokens == 4) {
      freq_data$continent  <- gsub("^([^_]+_[^_]+)_.*", "\\1", file_name)
      freq_data$class <- ifelse(freq_data$value == 0, "other", gsub("^.*?_[^_]+_([^_]+)_.*", "\\1", file_name))
    } 
    else if (num_tokens == 3) {
      freq_data$continent  <- gsub("^(.*?)_.*", "\\1", file_name)
      freq_data$class <- ifelse(freq_data$value == 0, "other", gsub("^.*?_(.*?)_.*", "\\1", file_name))
    }
    
    freq_data$year <- as.integer(gsub("^.*?_.*_(\\d{4}).*", "\\1", file_name))
    total_cells <- sum(freq_data$count)
    freq_data <- freq_data %>%
      group_by(continent, year) %>%
      mutate(percentage = round((count / total_cells) * 100, 2)) %>%
      ungroup()
    result_df <- bind_rows(result_df, freq_data)
  }
}

Udział procentowy obszarów pozostałych był nieprawdziwy, ze względu na potrójne wystąpienia klasy “pozostałe” dla każdego roku i każdego kontynentu. W celu poprawnego obliczenia udziału procentowego, ramkę danych pogrupowano oraz zsumowano udział procentowy do 100, względem pozostałych prawdziwych klas.

result_df_updated <- result_df %>%
  group_by(year, continent, class) %>%
  mutate(total_percentage = sum(percentage)) %>%
  group_by(year, continent) %>%
  mutate(percentage = ifelse(class == "other", 100 - sum(percentage[class != "other"]), percentage)) %>%
  ungroup() %>%
  select(-total_percentage) %>%
  distinct(year, continent, class, .keep_all = TRUE)

Rzeczywista liczba pikseli klasy “pozostałe” została obliczona na podstawie różnicy zsumowanej liczby pikseli wszystkich klas oraz trzech klas głównych (obszary leśne, zurbanizowane, rolnicze). W celu prezentacji statystyk dla całego świata, dane te pogrupowano według roku i kontynentu, dodając informacje o zsumowanej liczbie pikseli i udziale procentowym.

continents <- c("africa", "asia", "europe", "north_america", "oceania", "south_america")
files <- list.files(pattern = "_.*_\\d{4}.tif$", full.names = TRUE)
files <- files[seq(1, length(files), by = 3)]
cont_proc <- c()

for (file in files) {
  raster_data <- rast(file)
  freq_data <- as.data.frame(freq(raster_data))
  total_cells <- sum(freq_data$count)
  cont_proc <- c(cont_proc, total_cells)
}

result_df_updated <- result_df_updated %>%
  group_by(continent, year) %>%
  mutate(count = ifelse(class == "other", cont_proc[match(continent, c("africa", "asia", "europe", "north_america", "oceania", "south_america"))] - sum(count[class %in% c("agro", "forest", "urban")]), count))

world_aggregated <- result_df_updated %>%
  filter(continent != "world") %>% 
  group_by(year, class) %>%
  summarise(count = sum(count),
            percentage = round((count/sum(cont_proc)) * 100, 2),
            continent = "world")

result_df_updated <- bind_rows(result_df_updated, world_aggregated)

Taki zestaw danych pozwolił na obliczenie zmian powierzchni poszczególnych klas dla każdego kontynentu dla całego badanego okresu.

result_df_updated <- result_df_updated[1:10,]
result_df_updated
count continent year class percentage
241287787 africa 1992 other 62.84
56765277 africa 1992 agro 14.78
85626299 africa 1992 forest 22.30
312359 africa 1992 urban 0.08
198747887 asia 1992 other 45.07
97441833 asia 1992 agro 22.10
143880631 asia 1992 forest 32.63
880366 asia 1992 urban 0.20
23817815 europe 1992 other 21.18
47150640 europe 1992 agro 41.93


Wyniki

Globalne zmiany powierzchniowe

Obliczenie zmian udziału powierzchni może ukazać nie tylko trendy poszczególnych klas pokrycia terenu dla każdego kontynentu, ale także najważniejsze momenty w historii, które dany trend zaburzały lub odwracały.

Długotrwały, stabilny trend wzrostowy obszarów leśnych jest zauważalny jedynie w Afryce. Od kilku lat występuje również w Ameryce Południowej. Oceania jako jedyny kontynent odznacza się trendem spadkowym od roku 1992.

Zmiany udziału procentowego obszarów leśnych w latach 1992-2020

Największa zmiana powierzchni obszarów leśnych, sięgająca ponad 500 tys. km2, nastąpiła w Ameryce Południowej. Zmiana ta jest związana przede wszystkim z długotrwałą wycinką Puszczy Amazońskiej. Największym jednorocznym spadkiem powierzchni charakteryzuje się rok 2004, w którym nastąpiło zmniejszenie obszarów leśnych w Ameryce Południowej o ponad 120 tys. km2. Zmiany powierzchni obszarów leśnych dla pozostałych kontynentów są znacznie mniejsze i w ciągu 28 badanych lat zmieniały swoje trendy.

Zmiany powierzchni obszarów leśnych w latach 1992-2020


Globalne zmiany złożoności przestrzennej

Obliczone metryki posłużyły do ukazania stanu złożoności przestrzennej, jaka występowała na całym świecie w każdym z badanych lat. Wizualnie zmiany te są zauważalne głównie w dużych odstępach czasowych, dlatego poniżej przedstawiono wartości metryk jedynie w dwóch skrajnych latach: 1992 i 2020 roku.

Dane w obiekcie metrics nie posiadały przypisanego układu współrzędnych, zatem w celu ich poprawnego wyświetlenia należało przypisać każdemu kontynentowi jego indywidualne odwzorowanie. Następnym etapem była rasteryzacja danych ze względu na niepoprawne wyświetlanie obiektów w formie poligonowej oraz zmiana układu współrzędnych na wspólny, docelowy układ (tutaj: WGS84).

metrics_map <- function(class, year, metric){
  layers_crs <- c("+proj=aeqd +lat_0=8.5 +lon_0=21.5 +x_0=5621452.01998 +y_0=5990638.42298 +datum=WGS84 +units=m +no_defs", # africa
                    "+proj=aeqd +lat_0=47 +lon_0=94 +x_0=4340913.84808 +y_0=4812712.92347 +datum=WGS84 +units=m +no_defs", # asia
                    "+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +datum=WGS84 +units=m +no_defs", # europe
                    "+proj=aeqd +lat_0=52 +lon_0=-97.5 +x_0=8264722.17686 +y_0=4867518.35323 +datum=WGS84 +units=m +no_defs", # north america
                    "+proj=aeqd +lat_0=-19.5 +lon_0=131.5 +x_0=6988408.5356 +y_0=7654884.53733 +datum=WGS84 +units=m +no_defs", # oceania
                    "+proj=aeqd +lat_0=-14 +lon_0=-60.5 +x_0=7257179.23559 +y_0=5592024.44605 +datum=WGS84 +units=m +no_defs") # south america
    
    cont_layers = list("africa", "asia", "europe", "north_america", "oceania", "south_america")
    
    df <- c()
    for (i in 1:6){
      data <- df_sf[df_sf$continent == cont_layers[i],]
      data <- st_set_crs(data, layers_crs[i])
      df[[i]] <- data
    }
    
    rast_list <- list()
    for (j in 1:length(df)){
      rast <- df[[j]]
      empty_raster <- raster(rast , res = 30000)
      rast  <- terra::rasterize(rast, empty_raster, field = metric)
      rast  <- projectRaster(rast , crs = "+proj=igh +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs", res = 30000)
      origin(rast) <- 0
      rast_list[[j]] <- rast
    }
    merged <- do.call(merge, rast_list)
    assign(paste0(metric, '_', class, '_', year), merged, envir = .GlobalEnv)
}
metrics_map("forest", "2020", "ent")

projected <- projectRaster(ent_forest_1992, crs = 4326)


Entropia brzegowa

Rozkład przestrzenny wartości entropii brzegowej dla obu lat jest wizualnie bardzo podobny. Największa różnorodność krajobrazu leśnego występuje na większości obszarów Europy, a także w południowo-wschodniej Azji oraz na terenach tajgi w północnej Kanadzie. Rozkład wartości uległ zmianie - nastąpił lekki przyrost wszystkich obszarów z wartościami entropii powyżej 0, co oznacza zwiększenie zróżnicowania struktury krajobrazu.

Wartości entropii brzegowej w latach: a - 1992, b - 2020

Najbardziej widoczną zmianą średniej wartości entropii w ciągu 28 lat charakteryzuje się Ameryka Południowa., która zanotowała znaczący wzrost wartości metryki. Pozostałe kontynenty cechuje brak istotnej zmienności zróżnicowania przestrzennego krajobrazów. W ciągu ostatnich kilku lat wystąpił zauważalny trend wzrostowy entropii brzegowej dla każdego z kontynentów.

Zmiana średniej wartości entropii brzegowej w latach 1992-2020

Największym wzrostem entropii brzegowej między 1992 a 2020 rokiem odznacza się Amazonia. Znaczące spadki są zauważalne również na obszernych terenach leśnych we wschodniej Rosji, a także północnej Nikaragui i południowo-zachodniej Kanadzie. Obszary charakteryzujące się wzrostem entropii brzegowej są dwukrotnie liczniejsze od obszarów, dla których wartość metryki zmalała. Największe spadki zróżnicowania przestrzennego krajobrazu odnotowano na wielu obszarach Ameryki Południowej oraz Afryki Środkowej, a także w Birmie oraz stanie Minnesota w USA.

Zmiana wartości entropii brzegowej pomiędzy 1992 a 2020 rokiem

Wszystkie kontynenty odznaczają się znaczną przewagą obszarów ze wzrostem entropii brzegowej od obszarów, gdzie uległa ona spadkowi. Największa liczba lokalnych krajobrazów, dla których entropia brzegowa uległa znacznemu wzrostowi, występuje w Ameryce Południowej. Afryka odznacza się największym udziałem procentowym obszarów, które nie uległy zmianie, Europa natomiast - najmniejszym.

Liczba lokalnych krajobrazów w poszczególnych zakresach zmiany wartości entropii brzegowej


Względna informacja wzajemna

Przestrzenne zróżnicowanie wartości względnej informacji wzajemnej pomiędzy rokiem 1992 i 2020 jest dość widoczne. Dla obszaru całego świata nastąpił spadek wartości metryki, a więc zmniejszyła się jednolitość struktury przestrzennej lokalnych krajobrazów. Największy spadek wartości względnej informacji wzajemnej jest najbardziej widoczny na obszarze Amazonii, największy wzrost natomiast - w północnych Indiach.

Wartości względnej informacji wzajemnej w latach: a - 1992, b - 2020

Najbardziej widoczna zmianą średniej wartości względnej informacji wzajemnej wystąpiła w Ameryce Południowej w tym samym roku (2004), w którym jej średnie wartości entropii brzegowej znacznie wzrosły. Kontynentem znacznie odstającym od reszty jest Afryka, ze średnią wartością względnej informacji wzajemnej powyżej 0,6. W ciągu ostatnich kilku lat, w kontraście do rosnącej entropii, względna informacja wzajemna ulega zmniejszeniu na każdym z kontynentów.

Zmiana średniej wartości względnej informacji wzajemnej w latach 1992-2020

Największym spadkiem względnej informacji wzajemnej między 1992 a 2020 rokiem odznacza się Amazonia. Znaczące spadki wystąpiły również w Argentynie, środkowej części USA oraz na granicy Sudanu Południowego i Demokratycznej Republiki Konga. Łącznie na świecie więcej obszarów odznaczyło się spadkiem względnej informacji wzajemnej, aniżeli wzrostem. Największy wzrost obszarów, które w ciągu 28 lat stały się bardziej jednolite, nastąpił w północno-wschodnich Indiach, wschodnich Chinach oraz w południowo-zachodniej części Ameryki Południowej.

Zmiana wartości względnej informacji wzajemnej pomiędzy 1992 a 2020 rokiem

Europa jest jedynym kontynentem, którego większość obszarów odznacza się wzrostem względnej informacji wzajemnej. Pomimo dużego spadku jednolitości struktury przestrzennej na większości obszarów w Ameryce Południowej, wzrost względnej informacji wzajemnej również wystąpił na znacznej liczbie obszarów i jest relatywnie niewiele mniejszy.

Liczba lokalnych krajobrazów w poszczególnych zakresach zmiany wartości względnej informacji wzajemnej


Lokalne zmiany złożoności przestrzennej

Badanie zmian złożoności przestrzennej na poziomie lokalnym może znacznie przybliżyć przyczyny występujących zdarzeń i trendów, które wpłynęły na zwiększenie zróżnicowania struktury przestrzennej krajobrazów lub zwiększenia jej jednolitości.

Entropia brzegowa

Najwyższy spadek entropii brzegowej

Obszary z najwyższym spadkiem entropii brzegowej występują w południowej oraz wschodniej części Puszczy Amazońskiej w Brazylii. W ciągu 28 lat ponad 6 mln km2 obszarów przedstawionych na poniższej rycinie uległo zwiększonemu rozdrobnieniu krajobrazu. Tylko na ponad 600 tys km2 obszarach entropia brzegowa uległa lekkiemu spadkowi.

Obszary z najwyższym spadkiem entropii brzegowej - Brazylia



Największy wzrost entropii brzegowej

Znaczący wzrost jednolitości struktury krajobrazu wystąpił na granicy Sudanu Południowego oraz Demokratycznej Republiki Konga w Afryce. Zmiany te dotyczą prawie 300 tys. km2 obszarów, dla których nastąpił spadek entropii brzegowej - wystąpiło zalesienie wielu obszarów, które przedstawiono poniżej. Obszar ten wyróżnia się na tle terenów sąsiadujących, na których entropia brzegowa uległa wzrostowi.

Obszary z najwyższym wzrostem entropii brzegowej - Sudan Południowy oraz Demokratyczna Republika Konga



Względna informacja wzajemna

Najwyższy spadek względnej informacji wzajemnej

Oprócz najwyższego spadku zróżnicowania struktury krajobrazu, najwyższy spadek jednolitości struktury również występuje w Puszczy Amazońskiej. Obszary te są natomiast rozproszone po całej Amazonii, bez znacznego skupiska w Brazylii. W ciągu 28 lat ponad 5,5 mln km2 zanotowało spadek względnej informacji wzajemnej, co związane jest z długoletnią wycinką postępującą na obszarze Amazonii. W sąsiedztwie występuje również ok. 2 mln km2 obszarów, na których nastąpił lekki wzrost metryki.

Obszary z najwyższym spadkiem względnej informacji wzajemnej - Puszcza Amazońska



Największy wzrost względnej informacji wzajemnej

Obszary ze znaczącym wzrostem jednolitości struktury krajobrazu znajdują się na północy Indii i są związane z ubytkiem obszarów leśnych. Główne zmiany dotyczą ok. 150 tys. km2, na których zwiększyła się wartość względnej informacji wzajemnej.

Obszary z najwyższym wzrostem względnej informacji wzajemnej - Północne Indie

Identyfikacja obszarów o podobnej złożoności przestrzennej

Podsumowanie

W ciągu 28 badanych lat na zdecydowanej większości obszarów leśnych na świecie nastąpiła zmiana złożoności przestrzennej. Łącznie na świecie znajduje się znacznie więcej obszarów, na których różnorodność struktury krajobrazu uległa zwiększeniu, niż obszarów, gdzie wystąpił ich wzrost jednolitości. Świadczy to o postępującym rozdrabnianiu krajobrazu leśnego, co w ostatnich latach postępuje najmocniej. Głównym obszarem największych zmian w różnorodności oraz jednolitości struktury krajobrazu jest Puszcza Amazońska, gdzie od wielu lat występuje wycinka drzewostanu leśnego, która w ostatnich latach spowolniła i jest rekompensowana przez zalesianie.